%\documentclass[12pt,twoside]{report}
\documentclass[11pt]{report}
% note that the document can be single or double sided.
\usepackage{suthesis-3e}

\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{float}
\usepackage{xcolor}
\usepackage{listings}

\lstset{language=C++}

\lstset{breaklines}
\lstset
{
numbers=left,
framexleftmargin=10mm,
frame=none,
backgroundcolor=\color[RGB]{245,245,244},
keywordstyle=\bf\color{blue},
identifierstyle=\bf,
numberstyle=\color[RGB]{0,192,192},
commentstyle=\it\color[RGB]{0,96,96},
stringstyle=\rmfamily\slshape\color[RGB]{128,0,0},
showstringspaces=false
}


\setlength{\parindent}{0pt}




%\documentclass[11pt,a4paper,english,etoc,normalcite,sloppy]{zjuthesis}
%\definecolor{dgreen}{rgb}{0,.6,0}
%\newcommand\comments[1]{\textcolor{blue}{#1}}% my comments
%\renewcommand{\baselinestretch}{1.5}
%\input{preamble}


\includeonly{
time
%REPORT,
%Introduction,%
%%Background,%
%Metholody,
%Future_work,%
%body/Experiments,%
%body/Conclusion
 }


\DeclareMathOperator*{\argmax}{argmax}

\begin{document}
%type of this report
\reporttype{Early Stage Review Report}

%which department you are reporting to
\department2report{Department of Chemical Engineering}

%the name of your unversity
\reportUniversity{Imperial College London}

\title{Multi-scale simulation of multiphase multi-component flow in porous media using Lattice Boltzmann simulations}
%\ctitle{Xiaoyong WEI}
\author{Jianhui YANG \\(00559375)}

%the name of your department
\dept{Department of Chemical Engineering}

%the name of your supervisor
\principaladviser{Edo Boek}
\coprincipaladviser{Geoffrey Maitland}
%\submitdate{Jun 2 2011}



\pagenumbering{arial}

    \beforepreface
\input{abstract}
    \afterpreface



%\input{acknowledgement}

\chapter{Introduction}\label{chap1}
The rapidly increasing demand for fossil fuel pushes the petroleum industry to get a deeper understanding of earth systems. To better understand the earth system, we are interested in determining and predicting the transport properties within the basic elements of earth systems, the rocks. In order to determine, or predict the transport properties within the rocks, the underlying microscopic physics should be properly understood. Therefore we are keen to understand, simulate and predict pore-scale physics in porous media. However, this is normally difficult. Firstly, the geometry of the rocks can be highly heterogeneous and complicated. What's more, the transport phenomena inside are normally nonlinear and coupled. As an example, we consider the flow of water and oil in the rocks. 
The transport properties of some porous media have been determined by experiments in laboratories. However, they are time consuming and very costly. For example, the experimental measurement of relative permeabilities takes a long time (several months for a single rock core) and is therefore very expensive (costs are of the order of several 10$^5$ USD) in commercial petrophysics labs. The quality of measured data depends on the experience of operators. Some empirical relations were proposed to estimate the transport properties within the rocks \cite{walsh_1984,mavko_1998}. Due to the complexity of the structure, these empirical relations are not always effective, especially for multi-phase and multi-component flows. With the rapid development of the computer, numerical simulation is rapidly becoming an alternative solution for this problem. Several commercial enterprises already offer ``digital rock physics'' services as a complement or even alternative to experimental permeability measurements \cite{www_ingrain}. However, this problem is far from solved: the complex geometry of some classes of porous media, such as carbonate rocks, makes it difficult to model and predict the transport properties. Several numerical methods have been proposed to determine the transport properties of porous media, including Finite Element Method (FEM)\cite{kandhai_1998,saeger_1995}, Finite Differential Method (FDM)\cite{ladd_1994,adler_1990} and Network Modeling\cite{fatt_1956,dong_2009}. 
Network Modeling is an attractive option as it uses a simplification of the actual pore space geometry and therefore provides relatively simple numerical solutions \cite{blunt_2001,blunt_1991}. Network models work well for relatively simple rocks such as sand stones. However, for more complex porous media, such as carbonates, which have a very broad and multi-modal pore size distribution, the generation of network models is generally difficult \cite{knackstedt_2007}. For this reason, calculation of the transport properties directly on pore space images obtained from x-ray micro tomography (XMT)  has recently seen a tremendous development.
The direct methods mentioned above, FEM and FDM, generally require extensive meshing to obtain reasonable accuracy. They often are difficult to implement and are of low efficiency. Moreover, multi-physics processes including multi-phase and multi-component flow in complex geometries are very challenging for these traditional CFD methods.  As a powerful alternative, the lattice Boltzmann method (LBM) has been introduced as a novel CFD technique which is able to handle extremely complex geometries without simplification \cite{cancelliere_1990}. The Navier-Stokes equations with complicated boundaries can be solved accurately by the LB  method. These advantages match our requirement for the determination of transport properties in the pore space of the rocks. The LB method is easy to implement compared with other CFD algorithms. Most of the computation is local and thus it is ideal for parallel implementation. Due to the statistical physics background, it is easy to simulate multi-physics processes including multi-phase and multi-component flow, evaporation, condensation as well as cavitation. These advantages make the LB method an ideal numerical tool to study flow properties in porous media. %\\
The aim of this PhD study is to develop a lattice Boltzmann Method simulator for the transport phenomena associated with the flow of reservoir fluids, comprising of concentrated brines, $CO_2$ and hydrocarbons at the pore scale. It should be able to determine the transport properties in realistic rock geometries. This simulator should be reliable, efficient and accurate compared with other existing numerical methods and more importantly have the capability of handling complex geometries. We will use upscaling algorithms to predict flow properties at the core scale. At the end of the project, we aim to have a working LB simulation code to help design $CO_2$ re-injection operations in Qatari oil and gas fields. I would like to acknowledge that this PhD project is funded by Qatar Carbonates and Carbon Storage Research Centre (QCCSRC).
  

\chapter{Background}
\section{Literature Review}
Porous media consist of pore space and a solid matrix. The pores are generally interconnected which allows the fluids to flow and mass transfer to take place. The transport properties within a porous medium can be determined accurately if the physical equations can be solved. However, since the geometry is often extremely complicated and the problems in the porous media are nonlinear and coupled \cite{Keehm_PHD}, it is not practical to achieve an analytical solution. The porous medium is normally simplified as a homogeneous medium, and the parameters of the transport properties are measured by experiments in laboratories at the macroscopic scale. The experiments are very time consuming and costly. As computer power develops, several numerical methods have been proposed as an alternative solution for the estimation of transport properties within a  porous medium.\\

Network modeling was first proposed by Fatt \cite{fatt_1956} to investigate the capillary pressure. This method has been widely used to estimate the flow and transport in porous media in the past 20 years. It simplifies the pore space into nodes that represent pore junctions and bonds that connect pore nodes. This simplification keeps a certain amount of the geometry information at the pore scale including the connectivity and pore size distributions. At the same time it reduces the complexity of the geometry. A limitation of the network modelling method is that the algorithm that separates  the pore space into pore junctions and bonds is not unique. The distinction between the pores and the throats is not always clear in  reality. This makes the procedure of dividing somewhat arbitrary. The simplified network model cannot always reflect the real geometry, and therefore this affects the accuracy of the estimation. \\

More precise network models have been developed in the past 15 years to improve the performance. Maximal ball algorithms were proposed to construct the largest spheres centered on each void voxel that just fits in the pore space \cite{silin_2006}. 
%Okabe and Blunt (2004) used multiple point statistics to reconstruct 3D images using 2D image data and  successfully described the long-range connectivity of both sand stone and carbonate samples. 
A network model extracting process using the Maximal ball algorithm was developed based on a 3D representation of the pore space, obtained by micro-CT scanning tools. Flow and transport can be computed within this idealised geometry \cite{dong_2009}.\\

Several conventional CFD techniques have also been used to study the transport properties of porous media. Adler et al.\cite{adler_1990} used an alternating-direction-implicit (ADI) finite differential method to study Newtonian flow in 3D Fontainebleau sandstones with different porosities. Comparison between simulated permeability and experimental data was carried out; although the general shape of the experimental permeability versus porosity curve was predicted in quite an accurate way, the calculated permeability differs by, at most, a factor of 5 from the measured one. A regular mesh can be used in the finite difference method, but accurate computation with complicated boundaries requires an intensive mesh. Saeger et al. \cite{saeger_1995} used the finite element method (FEM) to solve the Stokes equations in periodic porous media. The calculated permeability achieved reasonable agreement with computed permeability given by other authors where the relative difference was from 0.31\% to 15\% \cite{larson_1988}. However, a very complicated irregular mesh needs to be used in finite element simulations for porous media which limits the use of this method. \\


The lattice Boltzmann method \cite{doolen_1990,chen_1992} uses a different approach to simulate the fluid flow. The LB method describes the fluid system by interactions of fictitious particle groups which reside on the lattice nodes. These particle groups are much bigger than the real fluid molecules, but show the same behaviour in density and velocity as the real fluid at the macroscopic scales. A number of studies showed that the LB Method can recover the Navier Stokes equations at the macroscopic scale \cite{frisch_1986,ladd_1994,chen_1992b}. The LB method uses a unique way to simulate the flow. While the conventional CFD method discretizes the governing equation in a top-down approach, the lattice Boltzmann method recovers the governing equations from the rules for discretized models in a bottom-up approach. \\

Since all the operations are taking place on the lattice nodes in the LB method, the lattice model is essential. Several lattice models were proposed for the LB method \cite{chen_1992b,humieres_1986,higuera_1989}. D'Humieres et al. \cite{humieres_1986} and Chen et al. \cite{chen_1992b} showed that in order to have isotropic relaxation of the stress tensor, the grid should be formed in 4D and then be projected into 2D and 3D. The grid in 2D and 3D used in this work was proposed by Qian et al. \cite{qian1992}. The 3D lattice model contains 19 velocities with the name D3Q19; the 2D lattice model contains 9 velocities with the name D2Q9. They are the most used lattice models in lattice Boltzmann simulations.Two limits of the lattice Boltzmann method were proposed by Qian et al. \cite{qian_1997} and Ladd \cite{ladd_1994}. They showed that the lattice Boltzmann method can recover the Navier-Stokes equations for low Mach number $Ma$ and low Knudsen number  $Kn$. This means the grid spacing has to be much smaller than the characteristic length ($Kn << 1$) and the fluid velocity should be much lower than the speed of sound ($Ma<<1$). \\

In its simplest incarnation, the LB method uses particle distribution functions which relax back to equilibrium using a single relaxation time (SRT) model. This is known as the BGK approximation and will be detailed in chapter 3. To overcome the disadvantages of the BGK model, such as numerical instabilities,  
Multiple-relaxation-time (MRT) lattice Boltzmann was developed by d'Humieres \cite{dhumieres1992}. The main idea of MRT is using different relaxation time parameters for different moments of macroscopic quantities. The collision step is carried out in the moment space which involves density, energy, momentum energy flux, diagonal and off diagonal components of stress tensor. The streaming step, in which the particles propagate to the neighboring nodes, is still done in velocity space as in the common LB method. The MRT model can overcome several common defects of the LB method such as fixed Prandtl number \footnote{The Prandtl number $Pr$ is a dimensionless number which describes the ratio of viscous diffusivity to thermal diffusivity. It is defined as $Pr=\frac{C_p \nu}{k}$, where $C_p$ is the specific heat, $\nu$ is the dynamic viscosity, $k$ is the thermal conductivity} ($Pr$=1 for common LBM) and fixed ratio between the kinematic and bulk viscosities. The stability improvement by using the MRT scheme would reduce the computational effort by at least one order of magnitude while maintaining the accuracy of the simulations \cite{dhumieres2001}. The numerical stability of MRT has been studied by (Lallemand and Luo 2000)\cite{lallemand2000} in detail. Their analysis showed that the MRT scheme led to a large improvement of numerical stability compared with the single relaxation time LB model. \\


There are three multi-component LB models for multi-phase and multi-component fluid simulations. The first multi-component LB model was proposed by Gunstensen et al (1991)\cite{color1991} with the name colour gradient model. Two components represent two types of fluid with their own distribution functions, and follow their own evolution equation. They are named red particles and blue particles. The collision step includes self-interactions and cross-interactions with other types of particles. A colour function gradient was introduced to calculate the surface tension between different phases. To segregate the phases, mixing near the interface should be minimized. A procedure called recolouring is proposed for this minimising process. This procedure is a very time consuming step, and this model also has some numerical stability problems for high density ratio and large surface tension. We found that this model is capable of simulating immiscible binary fluids with viscosity contrast but the same density. \\


The free energy model was developed by Swift et al. \cite{freeenergy}. This model includes thermodynamic equilibrium functions of phases, and a term describing the surface tension is added to the equilibrium distribution function. This allows the free energy model to specify the surface tension more easily than other multiphase multi-component models. It is also a fully thermodynamically consistent binary fluid lattice Boltzmann model. We found that the free energy model is able to simulate binary fluid systems with viscosity contrast but the same density. \\

A pseudo potential lattice Boltzmann model was developed by Shan and Chen \cite{shanchen}. The principal characteristic of this model is an interaction force between particles that is introduced to have a consistent treatment of the equation of state for a non-ideal gas. We found that the Shan-Chen model is able to simulate a binary fluid system with density contrast, but the viscosity of fluids must be the same.\\



A variety of studies on the transport properties of porous media using the LB method has been carried out. The absolute permeability of a reconstructed geometry from Fontainebleau sandstone was calculated using the lattice Boltzmann method by Jin \cite{jin_2004}. The porosity of the sandstone sample was $39.8\%$. Good agreement with the empirical formula was obtained. However, the permeability calculation was observed to be dependent on the viscosity. The viscosity dependence poses a severe problem for computing the permeability. Because the permeability is a characteristic of the physical properties of the porous media, it should only be related to the geometry of porous media. This viscosity dependence can be reduced by using the MRT scheme \cite{Pan_2006}. Pan et al.\cite{pan_2004} use a Shan-Chen multi-component lattice Boltzmann model to simulate multi-component flow in an idealised porous medium which comprises spheres with identical radius. The capillary pressure head was measured with saturation. Acceptable agreement with experimental data was obtained. The densities of the binary fluids were the same, and a small viscosity ratio of 1.8 was used in the simulation. This value of the viscosity ratio is generally too low for our systems containing oil, water and  supercritical $CO_2$ (the viscosity ratio of water and supercritical $CO_2$ is around 20). The efficiency of their code seems to be rather low, as the biggest geometry computed is a domain of $128^3$ voxels. Ramstad et al. \cite{ramstad_2009} used a color gradient lattice Boltzmann model to simulate an immiscible binary fluid system with equal viscosity in a Bentheimer sandstone where the geometry was taken from X-ray microtomography. Good agreement between the experimental data and simulated relative permeability of the wetting phase was obtained. The simulation under-predicts the relative permeability of the non-wetting phase at high wetting phase saturations. The viscosity and the density were kept the same in this study.  \\

Since 3D simulations normally require a large amount of computation, the efficiency of the algorithm becomes essential. A number of different lattice Boltzmann method implementation schemes were proposed to improve the efficiency of the algorithm. A parallel implementation can significantly improve the efficiency of the lattice Boltzmann method code. Several partition strategies \cite{skordos_1995,giovanni_1994,amati_1997,noble_1995} were proposed to divide the computation domains into slices or boxes. These decomposition strategies are only efficient when the workload is distributed homogeneously over the lattice \cite{grunau_1993}. Kandhai et al. \cite{kandhai_1998} proposed a new approach based on the Orthogonal Recursive Bisection (ORB) method. The ORB method can be used to generate approximately balanced decompositions by taking into account the workload on each lattice point. It is found to be $12\%$ to $60\%$ more efficient compared to conventional parallelization strategies. \\

A standard lattice Boltzmann method stores distribution functions for both fluid nodes and solid nodes. This wastes a large amount of memory and computing time, since no operation is carried out on solid nodes. A sparse data structure is required to reduce the memory and computing source usage. Implementation schemes with the name SHIFT were proposed to improve the efficiency of data storage in the computation \cite{ma_2010}. The information of solid nodes is no longer stored in these schemes, and the information regarding adjacent nodes is compressed to save memory. Numerical tests were carried out on four rock samples with porosities ranging from $10\%$ to $38\%$ to verify the performance of these schemes. The results showed that a reduction of $36\%$ to $82\%$ memory was obtained \cite{ma_2010}. \\

\section{Objectives} 
To calculate the transport phenomena associated with the flow of reservoir fluids, comprising of concentrated brines, CO2 and hydrocarbons in porous media, a parallel lattice Boltzmann code will be developed in C++ to study the flow and transport phenomena in porous media in two dimensions and three dimensions at the pore scale. We will consider the sparse storage strategy proposed by Mattila et al. \cite{mattila_2007} to reduce the usage of memory and CPU resources. Since the reservoir fluids systems we are going to simulate (water, oil and supercritical $CO_2$) normally have a density ratio ranging from 1 to 4 and a viscosity ratio ranging from 1 to 20, the multi-phase and multi-component lattice Boltzmann model should be capable of simulating an immiscible binary fluid system with similar properties. In order to select the most appropriate multi-phase and multi-component lattice Boltzmann model for the reservoir fluids in porous media, three multi-phase and multi-component lattice Boltzmann models including the Shan-Chen pseudo potential model, the Free Energy model and the Colour Gradient model are to be implemented. Numerical tests will be carried out to evaluate the performance of these models, and the most appropriate one will be selected for our project. 


This code package is expected to be able to:

\begin{itemize}
\item Simulate single-phase flow within complex geometries
\item Calculate the absolute permeability of porous media
\item Simulate multi-phase multi-component flow with viscosity and density contrast, surface tension and wettability in porous media
\item Calculate the relative permeability of porous media with binary immiscible fluids.
\item Compute the capillary pressure of binary immiscible fluids in porous media.
\end{itemize}

To obtain high accuracy and efficiency for the simulations of transport phenomena associated with the flow of reservoir fluids, we aim to overcome several difficulties of the lattice Boltzmann method:

\begin{itemize}
\item Improve numerical stability of the lattice Boltzmann method for low viscosity fluid flow simulations. 
\item Eliminate viscosity dependence of absolute permeability calculations.
\item Eliminate high spurious velocities in multi-phase and multi-component lattice Boltzmann method simulations. 
\item Reduce significantly the computing time in 3D lattice Boltzmann simulations. 
\item Reduce constraints on the viscosity ratio and density ratio of immiscible binary fluid systems that can be simulated by the lattice Boltzmann method.
\end{itemize} 



When the implementation of the code is complete,  we will first calculate the flow properties in 2D micro-model geometries comprised of capillary channels and junctions, where the corresponding microfluidic experiments will be carried out in a related separate PhD project (E. Chapman). Two examples of the micro-model geometries are shown in Figure \ref{Junction}. Our multi-phase lattice Boltzmann method code is to be used to simulate the flow and interaction of air and oil in these micro-model geometries. We aim to capture the capillary rising phenomena observed in the experiments.

\begin{figure}[H]
\begin{center}
\includegraphics[width=2in]{Junction.eps}
\end{center}
\caption{Example of geometry of micro-models}
\label{Junction}
\end{figure}

Secondly we aim to calculate flow properties in 3D pore space images of real reservoir rock, obtained (in a separate PhD project by PhD candidate S. Shah) from multi-scale imaging techniques including confocal microscopy and x-ray micro-tomography. With these images, a binary 3D matrix will be generated to describe the geometry of the rocks. These geometries will be used in our 3D single and multi-phase multi-component flow simulators. The absolute permeability will be computed by our single-phase lattice Boltzmann method simulator and compared with the experimental data. We will compute the absolute permeability of geometries generated from images with different resolution, to explore the relationship between simulated permeability and resolution of geometry.  We will calculate the multi-phase multi-component flow of brines, hydrocarbons and $CO_2$ in reservoir rocks. We aim to calculate relative permeabilities at the pore scale and these results will be compared with experimental data. We will use the absolute permeability and relative permeability from simulations at the pore scale combined with upscaling algorithms to predict flow properties at the core scale. A working LB simulation code will be developed at the end of the project to simulate the flow of oil, brines and supercritical $CO_2$ in the reservoirs at pore scale and core scale. These simulations are to help design $CO_2$ re-injection operations in Qatari oil and gas fields.
  




\chapter{Methodology, Viability and Progress to date}
Concepts of the lattice Boltzmann method will be introduced in this chapter; part of the progress to date will be presented to demonstrate the viability of the lattice Boltzmann method for the simulation of transport properties in porous media at the pore scale.

\section{Kinetic theory and Boltzmann equation}
The Boltzmann equation was introduced by Ludwig Boltzmann to describe the behaviour of a real gas using a statistical distribution of the gas particles. The gas is composed of interacting particles that can be described by classical mechanics. But because the number of particles is extremely large ($\scriptsize{\sim}10^{23}$), a statistical treatment was introduced to describe the averaged behaviour of the particles.  \\

A function $f(x,\xi,t)$ is introduced to describe the averaged distribution of particles. It represents the number density of particles at time $t$ and position $x$ with velocity $\xi$. We consider the gas in a control volume $dV=[x, x+dx]$. An external force $F=ma$ is applied to the system. At time $t$, the number of particles with velocity $[\xi, \xi+d\xi]$ is $dN=f(x,\xi,t)d\xi dV$. After a time $dt$, if there is no collision, the location of these particles will be $x'=x+\xi dt$, the velocity will be $\xi'=\xi+adt$. Therefore, we have the following equations:

\begin{eqnarray}
f(x+\xi dt,\xi+d\xi,t+dt)d\xi dxV-f(x,\xi,t)d\xi dx=0
\end{eqnarray} 

or

\begin{equation}
\frac{\partial f}{\partial t}+\xi \cdot \nabla_{x}f+ a \cdot \nabla_{\xi}f=0.
\end{equation}

Collisions may result in a change of the velocity of particles, so that the number of particles in the control volume $[x, x+dx] \times [\xi, \xi+d\xi]$ and $[x', x'+dx] \times [\xi', \xi'+d\xi]$ (which means particles are within $[x,x+dx]$, and the velocity of particles are within $[\xi, \xi+d\xi]$) will not be the same any more. A collision term 
$\Omega (f)$ is introduced to account for the change of distribution of particles. For an ideal gas, the collision term can be described as:

\begin{equation}
\Omega(f(\xi))=\int [f' f_1' -ff_1]B(\theta, |V|)d \theta d\epsilon d\xi_1
\end{equation} 

where $f(x,\xi,t)$ and  $f'(x,\xi',t)$ are the number densities before and after the collision, respectively. $V=\xi_1-\xi$ is the relative velocity between two particles, $\theta$ is the angle between $\xi_1-\xi$ and the the line linking the centres of two particles, $\epsilon$ is the projection angle of one particle over another, $B(\theta, |V|)$ is a non-negative function related to the interactions between the particles.\\

Let us introduce $\psi(\xi)$ as an arbitrary function of  $\xi$. If we set $\psi=1,\xi,|\xi^2|$, we can prove that:

\begin{equation}
\int \Omega(f)\psi(\xi)d \xi =0
\label{collision invariant}
\end{equation}

Any function which satisfies equation \ref{collision invariant} is called collision invariant \cite{bird1995molecular}. 

\section{H-theorem}
The H-theorem was introduced by Boltzmann in 1872 to describe the increase in the entropy of an ideal gas in an irreversible process. The H function is defined as:

\begin{equation}
H(t)=\overline{lnf}=\frac{\int flnf d\xi}{\int fd\xi}=\frac{1}{n}\int flnf d\xi
\end{equation}

Boltzmann proved that the H function is a monotonically decreasing function with time:
\begin{equation}
\frac{\partial H}{\partial t} \le 0.
\end{equation}
When the H function reaches a minimum value (Equation \ref{equi_H}), the system reaches equilibrium:
\begin{equation}
\frac{\partial H}{\partial t} =0
\label{equi_H}
\end{equation}


\section{Maxwell distribution}
Maxwell derived the probability equilibrium distribution for the speed of gas particles in 1860. It gives us the probability of a particle's speed being near a specified value as a function of temperature, velocity, and mass.

\begin{equation}
f=n\frac{1}{(2\pi R_gT)^{2/3}}exp[-\frac{(\xi-u)^2}{2R_gT}],
\label{maxwell_distribution}
\end{equation}

where $R_g$ is gas constant and $T$ is the thermodynamic temperature.\\


\section{Boltzmann-BGK equation}
In 1954, Bhatnagar, Gross and Krook \cite{BGK1954} introduced an approximate expression $\Omega_f$ for the collision term $\Omega(f)$. They proved that a simplified expression $\Omega_f$ which replaces the collision term $\Omega(f)$ should satisfy the following two properties:

\begin{itemize}
\item For collision invariant $\psi=(m,m\xi,\frac{1}{2}m\xi^2)$, the equation
      	\begin{equation}
	\int \psi \Omega_f d\xi=0
	\label{BGK1}
	\end{equation}
       should be satisfied.
\item It should satisfy Boltzmann's H-theorem:
	\begin{equation}
	\int(1+lnf)\Omega_fd\xi \le 0
	\label{BGK2}
	\end{equation}
\end{itemize}

They obtained the simplified term with the idea that the collision will lead the system to its equilibrium distribution $f^{eq}$. The rate of change is proportional to the difference of $f^{eq}$ and $f$. The scale factor $\nu$ is a constant: 
%The approximation term $\Omega_f$ is obtained:

\begin{equation} 
\Omega_f=\nu[f^{eq}(x,\xi)-f(x,\xi,t)]
\end{equation}

The Boltzmann equation is then simplified to:
\begin{equation}
\frac{\partial f}{\partial t}+\xi \cdot \nabla_{x}f+ a \cdot \nabla_{\xi}f=\nu (f^{eq}-f)
\label{Boltzmann-BGK}
\end{equation}

Equation \ref{Boltzmann-BGK} is called the Boltzmann-BGK equation. Bhatnagar, Gross and Krook also proved that $\Omega_f$ satisfies the properties \ref{BGK1} and \ref{BGK2}. The equation includes a collision time $\tau=\frac{1}{\nu}$, also called relaxation time, indicating the time interval between two collisions.  The Boltzmann-BGK equation is given as:

\begin{equation}
\frac{\partial f}{\partial t}+\xi \cdot \nabla_{x}f+ a \cdot \nabla_{\xi}f=-\frac{1}{\tau}(f-f^{eq})
\end{equation}

The macroscopic fluid density, velocity and energy can be calculated from the microscopic distribution function:
\begin{eqnarray}
&\rho(x,t)=mn(x,t)=\int f(x,\xi,t)d\xi & \\
&n\mathbf{u}(x,t)=\int \xi (x,\xi,t) d\xi&\\
&nR_gT(x,t)=\frac{1}{D} \int (v-\mathbf{u})^2f(x,\xi,t)d\xi&
\end{eqnarray}

where $m$ is the particle mass, $D$ is the dimension of the space, $\rho$ is the macroscopic density and $T$ is the temperature. By applying a Chapman-Enskog expansion, the macroscopic equations for mass, momentum and energy can be derived from the Boltzmann equation. The bulk viscosity is derived as:
\begin{equation}
\nu=\frac{\tau R_g T}{m}
\end{equation}

The equation of state (EOS) relating pressure and density is given by:
\begin{equation}
p=nR_gT
\end{equation}


\section{Single-phase lattice Boltzmann method}
The lattice Boltzmann Method (LBM) is a special discretization of the Boltzmann-BGK equation. Discretization of space, velocity and time are carried out in LBM. This procedure greatly simplifies the original Boltzmann equation. The location of particle distribution functions (PDFs) in space is restrained on the nodes of the lattice grid, and the particle velocity is simplified into a very limited number of lattice velocities. The number of discrete velocities is not unique, but the LB model using this discrete velocity model should be able to recover the macroscopic equations (in our study the Navier-Stokes equations). We take a 2-D model as an example. This model which is proposed by Qian et al\cite{qian1992} is well known and widely used. It contains 9 velocities and is known as D2Q9. In LBM, we assume that all the particle distribution functions (PDFs) have the same uniform mass (normally taken as 1 for simplicity). The lattice unit ($lu$) and time steps ($ts$) are important length and time units in LBM. We only discuss uniform mesh in this chapter ($\Delta x= \Delta y$).


\begin{figure}[H]
\begin{center}
\includegraphics[width=4in]{d2q9.eps}
\end{center}
\caption{D2Q9 lattice and velocities}\label{d2q9}
\end{figure}

Figure \ref{d2q9} shows the discretized velocity space $\{ \mathbf{e}_i \} ,\quad (i=0..8)$. The lattice velocity can be written as:

\begin{equation}
\mathbf{e}=e\begin{bmatrix}
0 & 1 & 0 &-1 &0 &1& -1& -1& 1 \\
0 &0 &1 &0 &-1& 1& 1& -1& -1 
\end{bmatrix}
\end{equation}

where $e=\Delta x / \Delta y$ is the local lattice speed with a unit of $lu \cdot ts^{-1}$ (where $lu$ is the lattice unit for length, $ts$ is the time unit in LBM). The relation with the local speed of sound is given as $c_s=\frac{e}{\sqrt{3}}$.\\

The continuous distribution functions associated with velocity are written as $f_i(x,t),\quad (i=0..8)$. We can obtain the Lattice Boltzmann equations for the D2Q9 model (Single Relaxation Time BGK) as:

\begin{equation}
f_i(\mathbf{x}+\mathbf{e}_i \Delta t,t+\Delta t)=f_i(\mathbf{x},t)-\frac{f_i(\mathbf{x},t)-f^{eq}_i(\mathbf{x},t)}{\tau}
\end{equation} 

Collision of the particles can be considered as a relaxation process towards equilibrium. The equilibrium distribution function for the D2Q9 model is a truncated Maxwell-Boltzmann distribution (Equation \ref{maxwell_distribution}) and is defined as \cite{qian1992}:

\begin{equation}
f^{eq}_i(x)=w_i\rho (x)[1+3\frac{\mathbf{e}_i \cdot \mathbf{u}}{e^2}+\frac{9(\mathbf{e}_i \cdot \mathbf{u})^2}{2e^4}-\frac{3\mathbf{u}^2}{c^2}]
\end{equation}

where the weight coefficients for the D2Q9 model are:
\begin{equation}
w_i=\begin{cases}
4/9 & i=0 \\
1/9 & i=1..4 \\
1/36 & i=5..8
\end{cases}
\label{wi_value}
\end{equation}

The macroscopic transport equations for mass, momentum and energy can be derived from the Boltzmann equation using a Chapman-Enskog expansion \cite{ekexpansion}. The kinematic viscosity $\nu$ in the D2Q9 model is obtained as:
\begin{equation}
\nu=c^2_s(\tau-\frac{1}{2})\Delta t
\end{equation}

%Its unit are $lu^2ts^{-1}$. 
Note that $\tau>1/2$  for positive viscosity. Numerical difficulties can arise as $\tau$ approaches $1/2$. The pressure is given by the equation of state for an ideal gas \cite{LBMODELING}:
\begin{equation}
P=\frac{nRT}{V}
\end{equation}

In single phase LBM, $RT=c_s^2=\frac{1}{3}$, so that the pressure can be computed as:
\begin{equation}
p=\rho c^2_s
\end{equation}


To implement a Lattice Boltzmann simulation, four major steps should be included in the code:
\begin{itemize}
\item Initialisation of distribution function $f_i(\mathbf{x},0)$
\item Collision step
	\begin{equation}
	f'_i(\mathbf{x},t)=f_i(\mathbf{x},t)-\frac{f_i(x,t)-f^{eq}_i(x,t)}{\tau}
	\label{collision_law}
	\end{equation}
\item Streaming step
	\begin{equation}
	f_i(\mathbf{x}+\mathbf{e}_i\Delta t,t+\Delta t)=f'_i(\mathbf{x},t)
	\label{streaming_step_equ}
	\end{equation}
\item Computation of macroscopic hydrodynamic quantities
	\begin{eqnarray}
	&\rho(\mathbf{x},t)=\sum_{i} f_i(\mathbf{x},t)& \label{macro_cal_1}\\
	&\rho \mathbf{u}(\mathbf{x},t)=\sum_{i} \mathbf{e}_i f_i(\mathbf{x},t)&\label{macro_cal_2}
	\end{eqnarray}
\end{itemize}

where $f'_i(x,t)$ represents the value of the distribution function $f_i(x,t)$ after collision.

\section{Bounceback Boundary Conditions}

Bounceback boundary conditions play a major role in the LBM simulation due to their simplicity, versatility and powerful capability of dealing with extremely complex boundaries. This boundary condition is usually used at fluid-solid interfaces due to its correspondence to the no-slip condition. This boundary condition is illustrated in Figure \ref{bounceback}. The densities moving toward the solid are bounced back into the fluid domain along reversed incoming directions. In the D2Q9 model, as shown in Figure (\ref{d2q9}) the bounceback condition can be described in terms of equations as:

\begin{eqnarray}
f_2=f_4,&f_5=f_7, f_1=f_3, &f_6=f_8
\label{bouceback_rule}
\end{eqnarray}

The standard bounceback condition places the boundary on the lattice nodes. Although mass and momentum are conserved, the accuracy is first order, while LBM is of second order \cite{ekexpansion}. Inamuro \cite{inamuro1995} found that the error produced by single relaxation time LBM with bounceback condition is sufficiently small and of second order if the relaxation parameter $\tau$ is close enough to 2. The bounce back conditions can be used without any influence on the order of the LBM, if $\tau$ is chosen in the range (0.5,2). Bounceback conditions in multi-relaxation time LBM are always of second order and not dependent on the selection of relaxation parameters \cite{Pan_2006}. Furthermore, the bounce back condition is the most efficient one for arbitrarily complex geometries \cite{breuer2000}. Many researchers contributed to this ongoing discussion \cite{noble_1995,ziegler1993}. A second order scheme with the name Half-Way bounce back condition was proposed by Ziegler\cite{ziegler1993}. In this boundary condition, the surface is a solid boundary placed between two neighbouring lattice sites with the same distance $\Delta x /2$. It is illustrated in Figure (\ref{bounceback}). This half-way boundary condition has been implemented in our simulator. \\ 

The boundary condition can be integrated in the collision step. The particles in the fluid domain will follow the collision law (Equation \ref{collision_law}), whereas the particles in the solid domain follow the bounceback rule defined in Equation (\ref{bouceback_rule}). 


\begin{figure}[H]
\begin{center}
\input{bounceback.pstex_t}
\end{center}
\caption{Bounceback condition}\label{bounceback}
\end{figure}


\section{Periodic Boundary Condition}
The Periodic boundary condition could be the simplest boundary condition for LBM. In the periodic system, the fluid flow out through one face will reenter into the opposite face of the domain. So the edges of the simulation domain could be treated as if they are attached to the other side of the domain. For the boundary nodes, their neighbouring nodes are located at the opposite side of the boundary. (Fig \ref{pbc})

\begin{figure}[H]
\begin{center}
\input{pbc.pstex_t}
\end{center}
\caption{Periodic boundary condition}\label{pbc}
\end{figure}

Assuming that the length of computational domain is $L$, the periodic boundary condition could be described by the formula (taking the x direction as an example, the derivation for the y direction is straightforward):
\begin{eqnarray}
&f_i(0,y,z,t+\Delta t)=f'_i(L,y,z,t)&\\
&f_i(L,y,z,t+\Delta t)=f'_i(0,y,z,t)&
\end{eqnarray}


\section{Fixed Pressure or Velocity Boundary}
The fixed pressure or velocity boundary can be achived by updating the distribution function after the streaming step with an equilibrium distribution function that has the desired pressure/velocity. 

Guo et al \cite{guo_2002} proposed another approach to apply a pressure/velocity boundary. The distribution function after collision was decomposed into two parts; equilibrium and non-equilibrium: $f(x,t)=f^{eq}(x,t)+f^{neq}(x,t)$. The distribution functions of boundary points are updated using the non-equilibrium distribution from the neighbouring points and the equilibrium distribution value computed using the desired pressure/velocity. He showed that this scheme is numerically stable and of second order. 

\section{Multi-Relaxation-Time (MRT) scheme for the lattice Boltzmann method}
MRT allows independent adjustment of bulk
\footnote{Bulk viscosity, also called volume viscosity is important for simulations where fluid compressibility is essential. It appears in the compressible Navier-Stokes equation: 
\begin{equation}
\rho(\frac{\partial v}{\partial t}+v \cdot \nabla v)=-\nabla \rho +\mu\nabla^2 v+f+\mu^v\nabla(\nabla \cdot v)
\end{equation}
where $\mu^v$ is the bulk viscosity. In the incompressible Navier-Stokes equation, this term disappears because the divergence of the velocity of an incompressible fluid, $\nabla \cdot v$, equals 0} 
and shear viscosities which significantly improves the numerical stability for a low viscosity fluid. In the single relaxation time LB model, the collision term is relaxed by a single parameter $\tau$, while it could be relaxed instead by a matrix $\Lambda$:

\begin{equation}
f_i(x+c_i \Delta t)-f_i(x,t)=-\Lambda_{ij}[f_j-f_j^{eq}], \quad i=1,2,\ldots,b
\end{equation}

The matrix $\Lambda$ is a full matrix of constants. For the D2Q9 model, the transformation matrix $M$ is given as \cite{lallemand2000}:
\begin{equation}
M=\begin{pmatrix}
1  &1  &1  &1  &1  &1  &1  &1  &1  \\
-4 &-1 &-1 &-1 &-1 &2  &2  &2  &2\\
4  &-2 &-2 &-2 &-2 &1  &1  &1  &1\\
0  &1  &0  &-1 &0  &1  &-1 &-1 &1\\
0  &-2 &0  &2  &0  &1  &-1 &-1 &1 \\
0  &0  &1  &0  &-1 &1  &1  &-1 &-1\\
0  &0  &-2 &0  &2  &1  &1  &-1 &-1\\
0  &1  &-1 &1  &-1 &0  &0  &0  &0\\
0  &0  &0  &0  &0  &1  &-1 &1 &-1
\end{pmatrix}
\end{equation}

The single relaxation time LB model can be obtained by specifying $\Lambda$ as a diagonal matrix with identical value:

\begin{equation}
\Lambda_{ij}=\frac{1}{\tau}\delta_{ij}
\end{equation}

The macroscopic quantities are calculated in the same way as in the LBGK model. Instead of considering distribution functions, MRT employs several moments corresponding to macroscopic quantities and their flux. These quantities can be relaxed with different time scales. A matrix $M$ transforms the distribution functions $f_i$ from the distribution space to the moment space: 

\begin{equation}
m=M\cdot f, \quad f=M^{-1}\cdot m
\end{equation}

The moment space for the D2Q9 model is:

\begin{equation}
m=(\rho,e,e^2,j_x,q_x,j_y,q_y,p_{xx},p_{xy})^{T}
\end{equation}
where $e$ is the energy, $j_x,j_y$ are the momenta in $x$ and $y$ directions, $q_x,q_y$ are energy fluxes and $p_{xx},p_{xy}$ is the stress tensor. The collision is carried out in the moment space by multiplying the transformation matrix $M$; the left and right hand side of equation \ref{collisionmrt} can be transformed into the moment space as:

\begin{equation}
f_i'(x,t)=f_i(x,t)-\Lambda_{ij}[f_j-f_j^{eq}]
\label{collisionmrt}
\end{equation}
\begin{equation}
\mathbf{m}'=\mathbf{m}-\mathbf{S}[\mathbf{m}-\mathbf{m}^{eq}]
\end{equation}

where $\mathbf{m}^{eq}=\mathbf{M}\mathbf{f}^{eq}$ is the equilibrium equation in moment space. $S=M \Lambda M^{-1}=diag(s_1,s_2,\ldots,s_b)$. The corresponding relaxation time for moment $m_i$ is $s^{-1}_i$. After the collision step, the moment $m'$ is transformed back into distribution function space by multiplying $M^{-1}$ for the streaming step $f_i(x+e_i\Delta t,t+\Delta t)=f_i'(x,t)$ which will be carried out in the same way as in the  single relaxation time LB model. 

The relaxation parameters and equilibrium functions of the moments are:



\begin{equation}
S=(0,s_e,s_{e^2},0,s_q,0,s_q,s_{\nu},s_{\nu})
\end{equation}
\begin{equation}
m^{eq}=\rho(1,-2+3u^2,\alpha+\beta u^2,u_x,-u_x,u_y,-u_y,u^2_x-u^2_y,u_xu_y)^{T}
\end{equation}

where $\alpha$ and $\beta$ are adjustable parameters and may be chosen as $\alpha=1,\beta=-3$ \cite{dhumieres1992}. The kinematic viscosity and volume viscosity are given by
\begin{eqnarray}
&\nu=c^2_s(\frac{1}{s_{\nu}}-\frac{1}{2})\Delta_t& \\
&\zeta=c^2_s(\frac{1}{s_{e}}-\frac{1}{2})\Delta_t&
\end{eqnarray}

where $c_s$ is the local sound speed that equals $1/\sqrt{3}$ in the single phase LB model, $s_{\nu}$ and $s_e$ are parameters given by users to set the kinematic viscosity and volume viscosity. 


\section{Verifications for the single-phase lattice Boltzmann method}

\subsection{Poiseuille flow simulation and flux calculation in a narrow channel}

Poiseuille flow in a channel with two parallel solid surfaces was studied. This model can be regarded as a simplified pore structure in porous media. Periodic boundaries are used in the flow direction. The analytical velocity profile in a slit of width $2a$ is parabolic and is given by the Poiseuille equation \cite{chin_2002}:

\begin{equation}
u(x)=\frac{G^*}{2\nu}(a^2-x^2)
\label{poisseuille_analytical}
\end{equation}

where $G^*$ is a pressure gradient or a body force applied on the fluid. The flux across the outlet can be computed by $\int u(x) dx$. We compute the flux across the outlet of slits with different width by using the Half-Way and Standard bounce back boundary conditions for the walls, while periodic boundary conditions are used for the inlet and the outlet boundaries. The width varies from 1  to 51 lattice points. The simulated and analytical velocity profiles are shown in Figure \ref{velocity_flux}. Excellent agreement is obtained in the simulations for all widths. Even when the slit has only one cell spacing for the fluid, the LBM code still gave excellent results for the velocity estimation. This excellent agreement is due to the use of the half-way bounce back boundary condition for the walls along the flow direction and the use of the Multiple-Relaxation-Time scheme (MRT). The accuracy of the half-way Bounceback boundary condition is second order\cite{LBMODELING}, and the Poiseuille flow is a second order problem (Analytical solution is shown in Equation \ref{poisseuille_analytical}). Therefore LBM can predict the velocity field perfectly. The MRT scheme eliminates the viscosity dependence of the velocity calculations. A standard lattice Boltzmann method simulation without the MRT scheme introduced an error in the velocity ranging from $8\%-17.8\%$ depending on the viscosity of the fluid\cite{Pan_2006}. In our results in Figure \ref{velocity_flux}, accurate velocity fields can be obtained with arbitrary viscosity values.\\ 

\begin{figure}[H]
\begin{center}
\includegraphics[width=7.2in]{velocity_flux.eps}
\end{center}
\caption{Simulated velocity compared with analytical solution for different channel widths of 19, 9, 4, 3, 2, 1 $lu$}\label{velocity_flux}
\end{figure}


The flux is calculated by integrating the velocity along $x$. The simulated flux is calculated by summing the velocity along the direction that is perpendicular to the flow direction. Since we know from our previous simulation that the velocity profile is accurate, the calculation of the flux becomes a procedure of numerical integration. A midpoint rule is used to calculate the flux. Therefore, the more integration points we have, the higher accuracy we can obtain. This analysis matches our numerial experiment (Fig.\ref{error_flux}). When the channel spacing is only one cell, although we can predict the velocity accurately, the error in the flux calculation is 50\%. However, when the channel is 3 lattice sites wide, the error in the flux calculation decreases dramatically to 5\%, which is an acceptable accuracy. The error will continue to decrease as the width of channel increases. We can obtain very good results with a width bigger than 30. In this numerical experiment, we also learn that the half-way bounce back boundary condition introduces a significantly lower error in flux estimation, compared with standard bounceback boundary conditions.% To calculate the permeability of a geometry with pore structures accurately, we need to keep the minimum pore diameter bigger than 3 lattice in order to include enough integration points for flux calculation. 

\begin{figure}[H]
\begin{center}
\includegraphics[width=4in]{tube_flux_error.eps}
\end{center}
\caption{Error of simulated flux of a slit}\label{error_flux}
\end{figure}


\subsection{Flow through A Pipe}
A simple 3D cylindrical tube model with a grid of 60x60x60 is studied with our code. Figure \ref{pipe_flux_geo} shows the geometry of the pipe model. The porosity is dependent on the radius. The fluid velocity has a parabolic profile similar to Poisseuille flow. The analytical solution for this laminar Poiseuille flow is\cite{Keehm_PHD}:

\begin{equation}
Q=\frac{\pi R^4}{8\eta}\frac{\partial P}{\partial x}
\end{equation}


\begin{figure}[H]
\begin{center}
\includegraphics[width=3in]{Pipe_flux_geometry.eps}
\end{center}
\caption{The Pipe model}\label{pipe_flux_geo}
\end{figure}




where $\eta$ is the dynamic viscosity, $R$ is the radius of the tube and $\frac{\partial P}{\partial x}$ is the pressure gradient. The permeability $\kappa$ is computed by applying Darcy's Law:

\begin{equation}
\kappa=\frac{\pi R^4}{8 A}
\end{equation}

where $A$ is the cross-sectional area. The calculated permeability from our code and Keehm's paper \cite{Keehm_PHD} are shown in Figure \ref{tube_flux}. Good agreement with the theoretical prediction is obtained by both our code and Keehm's results. Relative errors in the percentage of the analytical solution with increasing radius of the tube are shown in Figure (\ref{tube_per_error}). The results show the trend of error decreasing from $40\%$ to around $5\%$  with tube radius increasing. This is because a big radius gives a finer mesh on the boundaries. My code also gives a lower error of $2\%-5\%$ compared with Keehm's method. These better results are probably due to the use of the MRT scheme.\\ 

Another interesting result can be observed in Figure \ref{tube_per_error}: the error observed in a simulation with an integer radius is bigger than that for a non-integer radius. This reduction of error is due to the boundary condition. The input file of the boundary is a binary matrix with 0 and 1 that represents pore nodes or solid respectively. If the radius is one half larger than an integer value, the actual boundary will be located between the pore nodes and solid nodes. This introduces an effective half-way bounce back boundary condition. Otherwise, a standard bounce back boundary condition (SBB) is used. We have already seen that half-way bounce back boundary conditions have a second order accuracy which is one order higher than the standard bounce back boundary condition. This numerical experiment suggests that the use of half-way bounceback boundary conditions allows more accurate calculation of the permeability. We can also learn from these simulation results that boundary conditions affect the permeability calculation significantly\cite{sengupta2011error}. We will continue studying this topic in the future.

\begin{figure}[H]
\begin{center}
\includegraphics[width=4in]{tube_flux.eps}
\end{center}
\caption{Simulated permeability from our code, Keehm's paper and the theoretical prediction as a function of tube radius, both are in lattice unit}\label{tube_flux}
\end{figure}

\begin{figure}[H]
\begin{center}
\includegraphics[width=4in]{tube_per_error.eps}
\end{center}
\caption{Relative Error in percentage of the analytical solution with the radius of the tube, my results are shown in black star points, Keehm's results are shown in blue triangle points}\label{tube_per_error}
\end{figure}

\subsection{Simulation of flow in fibrous porous media}

A porous medium with rectangular periodic arrays of cylinders with elliptical cross section \cite{yang_2000} is investigated to evaluate the performance of our code. The square packing configuration for rows of elliptical fibers is shown in Figure \ref{fibrous_geo}. The LBM simulation domain is shown in Figure\ref{single_porous}. Periodic boundary conditions are applied in x,y and z directions. The porosity is determined by the parameters $Lx,Ly,a,b$. An analytical formulation of the permeability is derived by Phelan et al.\cite{phelan_1996}:

\begin{equation}
K=\frac{2L_y}{L_x}\frac{1}{\int_{-L_y/2}^{L_y/2} \frac{h^3}{3}}
\end{equation}

where $h$ is a function of the y-coordinate, as shown in Figure \ref{single_porous}. It indicates the distance between the cylinder and the computational boundary.To evaluate our simulated permeability with this theoretical prediction, a dimensionless permeability is defined as:

\begin{equation}
K^*=\frac{4K}{ab}
\end{equation}

\begin{figure}[H]
\begin{center}
\includegraphics[width=1.0in]{fibrous_porous_geo.eps}
\end{center}
\caption{Geometry of fibrous porous media}\label{fibrous_geo}
\end{figure}

\begin{figure}[H]
\begin{center}
\includegraphics[width=3in]{single_porous.eps}
\end{center}
\caption{LBM simulation domain}\label{single_porous}
\end{figure}

The permeability for a single cylinder with periodic boundary conditions is calculated using a 150x150 2D mesh in the simulation. The results are compared with Keehm's LB results \cite{Keehm_PHD} and the analytical solution. It can be seen from Figure \ref{elliptical_per} that good agreement with the analytical solution is obtained. For porous media  with a low porosity 
(less than $23\%$), the prediction of the permeability from both Keehm's  and our results showed a big error (up to $47\%$ for my results, $80\%$ for Keehm's results). This is due to the lack of fine mesh between adjacent cylinders. 

\begin{figure}[H]
\begin{center}
\includegraphics[width=4in]{elliptical_cylinder_permeability.eps}
\end{center}
\caption{Simulated dimensionless permeability of arrays of cylinders with elliptical cross sections}\label{elliptical_per}
\end{figure}



%\subsection{Permeability simulation on Clashach Sandstone and Doddington Sandstone}

%A Clashach sandstone with a porosity of $11.96\%$ and a Doddington sandstone with a porosity of $18.72\%$ are evaluated for permeability simulations.The aim of these simulations is to explore the performance of our LB code for real rock geometries. Because the measurements of permeability of rocks are very time consuming in the labs, a reliable simulator can significantly reduce this measuring time. \\


%The pore geometry of the rocks were obtained by Al-Ansi et al(2011) using X-ray microtomography with five different resolutions: $5.999\mu m/pixel,8.999\mu m/pixel,11.998\mu m/pixel,14.998\mu m/pixel,19.997\mu m/pixel$. A sample with resolution of 19.997 $\mu m/pixel$ is shown in Figure(\ref{Dodding20GEO}), where the red part represents solid,while the grey part represent pores. Carrying out simulations on geometries with different resolutions can help us investigate the digitization effect.\\


%\begin{figure}[!hpg]
%\begin{center}
%\includegraphics[width=3in]{Doddington20.eps}
%\end{center}
%\caption{Dodding Sandstone with a resolution of 19.997 $\mu m/pixel$}\label{Dodding20GEO}
%\end{figure}

%In Figure(\ref{Clashach}), we show the simulated permeability of Clashach Sandstone as a function of resolution. The measured permeability in the lab ,simulated permeability from lattice Boltzmann method, Navier-Stokes simulations and network simulations are showed as a solid line ,squares,triangles and stars respectively. In general, the simulations can predict the permeability well only for low resolution samples. Big errors were observed in high resolution simulations.The lattice Boltzmann method and network simulations seem to give better results than the Navier-Stokes simulation. \\

%The simulation of low resolution samples (19.997 $\mu m/pixel$) can predict the permeability with the highest accuracy. Lattice Boltzmann method and network simulation obtained good agreement with experimental data. However, I found that the simulated permeability increase with increasing resolution. All simulations including lattice Boltzmann method, Navier-Stokes and network seem to over predict the permeability of high resolution samples.

%\begin{figure}[!hpg]
%\begin{center}
%\includegraphics[width=3in]{DATA_Clashach.eps}
%\end{center}
%\caption{Simulated permeability and experimental data of Clashach Sandstone}\label{Clashach}
%\end{figure}


%We investigated the effect of the quality of mesh on the simulation results. Single phase flow was simulated with a geometry of Doddington Sandstone in different resolutions ($5.999\mu m/pixel$,$8.999\mu m/pixel$,$11.998\mu m/pixel$,$14.998\mu m/pixel$,$19.997\mu m/pixel$). A regular cubic mesh was used in the simulation. Simulations with the original mesh, the mesh that was refined once, and the mesh that was refined twice were carried out to calculate the absolute permeability. The results are shown in Figure(\ref{Dod}). We find that the permeability from high resolution geometries is more accurate than that in low resolutions. And the quality of the mesh affects the value of simulated permeability significantly. Mesh refinement can significantly improve the accuracy of the calculated permeability. 
%If the mesh is fine enough, the lattice Boltzmann method can predict the permeability with reasonable accuracy.We are currently trying to rationalise the resolution effects observed for both rock types.

 
%\begin{figure}[!hpg]
%\begin{center}
%\includegraphics[width=3in]{DATA_Doddington.eps}
%\end{center}
%\caption{Simulated permeability and experimental data of Dodding Sandstone}\label{Dod}
%\end{figure}

\section{Multi-component lattice Boltzmann method}
A multi-component system consists of separate
chemical components such as oil and water. Due to their economic
importance, such systems have been studied extensively. Three multi-component lattice Boltzmann models have been implemented, and their performance has been evaluated with some numerical experiments. 

\subsection{Shan-Chen pseudo potential model}
A pseudo-potential model for multi-phase and multi-component lattice Boltzmann simulation was proposed by Shan and Chen \cite{shanchen}. The principal characteristic of this model is an interaction force between particles which is introduced to have a consistent treatment of the equation of state for a non-ideal gas. For the D2Q9 model, the attractive force $F$ between the particles is described as:

\begin{equation}
F(\mathbf{x},t)=-G\Psi (\mathbf{x},t)\sum_{i=1}^{8}w_i \Psi (\mathbf{x}+\mathbf{e}_i\Delta t,t)\mathbf{e}_i
\label{Shanchenforce}
\end{equation}

where $G$ is a parameter determining the interaction strength between neighbouring particles, $w_i$ are weightings defined in Equation (\ref{wi_value}). It also determines whether the interaction is attractive or repulsive. In multi-component immiscible flow, same component particles will attract each other while particles from different components will repel each other. From this, phase separation can be obtained.  $\Psi$ is the interaction potential of the form:

\begin{equation}
\Psi (\rho)=\Psi_0exp(\frac{-\rho_0}{\rho})
\label{potentialfunction}
\end{equation}

where $\Psi_0$ and $\rho_0$ are constant parameters. Figure (\ref{potential})\cite{LBMODELING} shows the interaction potential function with $\Psi_0=4$ and $\rho_0=200$, where parameters can be selected arbitrarily. When $G$ is chosen as a negative value, the interaction between particles is attractive. The strength of the potential function depends on the density: the interaction strength increases with the density. Therefore the liquid, which has a higher density than the vapour, experiences a stronger cohesive force, which leads to surface tension phenomena \cite{LBMODELING}.

\begin{figure}[H]
\begin{center}
\includegraphics[width=3in]{potential_function.eps}
\end{center}
\caption{Interaction potential function of Shan-Chen model}\label{potential}
\end{figure}


\subsection{The Free Energy lattice Boltzmann model}
Swift et al. \cite{freeenergy} developed a new thermodynamically consistent
binary fluid LB model by introducing an equilibrium state associated
with a free energy functional, corresponding pressure tensors and
chemical potentials. A correct choice of collision rules ensures that
the system evolves toward minimization of the free energy
functional.

The thermodynamic properties of a binary fluid system can be
described by a Landau free energy functional \cite{freeenergy}:

\begin{equation}
\Psi_{FE}=\int_V (\psi_b+\frac{\kappa}{2}(\partial_{\alpha} \phi)^2)dv
+\int_S \psi_s ds
\end{equation}

where $\kappa$ is a parameter determining the strength of the surface tension $\sigma$ \cite{briant_2004}: $\sigma= \sqrt{8 \kappa A/9}$ , $\psi_b$ is the bulk free energy density and has the form:
\begin{equation}
\psi_b=\frac{c^2}{3} \rho ln\rho +A (-\frac{1}{2}\phi^2+
\frac{1}{4}\phi^4)
\end{equation}

$A=0.04$ is a parameter suggested by Pooley et al. \cite{pooley_2009}. $\phi$ is the order parameter representing the concentration of
components defined as:

\begin{equation}
\phi=\frac{\rho_a-\rho_b}{\rho_a+\rho_b}
\end{equation}

This free energy density function can lead to a phase separation
with $\phi=\pm 1$. The gradient term represents an energy
contribution and is related to surface tension. The second
integral is over the surface and describes the interactions between
fluids and the solid surfaces. $\psi_s$ is taken as:
\begin{equation}
\psi_s=-h \phi_s
\end{equation}

$\phi_s$ is the order parameter value at the solid surface.
Minimizing the free energy shows that the gradient in $\phi$ at the
solid surface is :

\begin{equation}
\kappa \partial_{\perp} \phi=\frac{d \psi_s}{d \phi_s}=-h
\label{wettingFE}
\end{equation}

The contact angle $\theta$ is related to the parameter $h$ and is
formulated as:
\begin{equation}
h=\sqrt{2\kappa A} sign(\frac{\pi}{2}-\theta)
\sqrt{cos(\frac{\alpha}{3})[1-cos(\frac{\alpha}{3})]}
\label{contact_angle_equation}
\end{equation}

where $\alpha=cos^{-1}(sin^2 \theta)$.

The pressure tensor $P_{\alpha\beta}$ and the chemical potential $\mu$ are defined in the usual way:

\begin{eqnarray}
&\mu=\frac{\partial \Psi}{\partial \phi}=A(-\phi+\phi^3)-\kappa
\nabla^2 \phi&\\
&P_{\alpha
\beta}=(P_b-\frac{\kappa}{2}(\partial_{\gamma}\phi)^2-\kappa \phi
\partial_{\gamma \gamma}\phi)\delta_{\alpha \beta}+
\kappa(\partial_{\alpha}\phi)(\partial_{\beta}\phi)& \\
&P_b=\frac{c^2}{3}\rho+A(-\frac{1}{2}\phi^2+\frac{3}{4}\phi^4)&
\end{eqnarray}

The hydrodynamics and thermodynamics of the binary fluids are
described by the Navier-Stokes and convection-diffusion
equations respectively:

\begin{equation}
\frac{\partial \phi}{\partial t}+\nabla (\phi u)= \mu \nabla ^2 \mu
\end{equation}


A new distribution function $g_i(x,t)$ is introduced to describe the
order parameter $\phi=\sum\limits_i g_i$ and is related to the
convection and diffusion. The distribution function $f_i(x,t)$ is
related to the fluid density and momentum as usual.\\

The time evolution equations using the MRT scheme can be described as:
\begin{eqnarray}
\text{Collision Step:} & f'_i(x,t)=f_i(x,t)-M^{-1}SM(f_i-f^{eq}_i) & \\
& g'_i(x,t) =g _i^{eq}(x,t)& \\
\text{Streaming Step:} & f_i (x+e_i \Delta t,t+\Delta t) =f'_i(x,t)&
\\
&g_i (x+e_i \Delta t,t+\Delta t) =g'_i(x,t)&
\end{eqnarray}

An appropriate choice of the equilibrium distribution functions can
reproduce the macroscopic equations in the continuum limit. In our case, we have used the equilibrium distribution functions proposed by Kusumaatmaja \cite{kusumaatmaja_phd}:

\begin{eqnarray*}
f_{i}^{eq}=&\frac{w_i}{c^2}(P_b-\kappa \phi\nabla^2 \phi+e_{i\alpha}\rho u_{\alpha}+\frac{3}{2c^2}[e_{i\alpha}e_{i\beta}-\frac{c^2}{3}\delta_{\alpha\beta}](\rho u_{\alpha}u_{\beta}+\lambda[u_{\alpha} \partial_{\beta}\rho+u_{\beta} \partial_{\alpha}\rho+\delta_{\alpha\beta}u_{\gamma}\partial_{\gamma}\rho])& \\
&+\frac{\kappa}{c^2}(w^{xx}_i \partial_x \phi \partial_x \phi+w^{yy}_i \partial_y \phi \partial_y \phi+w^{zz}_i \partial_z \phi \partial_z \phi+w^{xy}_i \partial_x \phi \partial_y \phi+w^{yz}_i \partial_y \phi \partial_z \phi+w^{zx}_i \partial_z \phi \partial_x \phi),&\\
&g_i^{eq}=\frac{w_i}{c^2}(\Gamma \mu+e_{i\alpha}\phi u_{\alpha}+\frac{3}{2c^2}[e_{i\alpha}e_{i\beta}-\frac{c^2}{3}\delta_{\alpha\beta}]\phi u_{\alpha u_{\beta}})&
\end{eqnarray*} 

where $\Gamma$ is a mobility parameter that sets the diffusivity of the interface and $w$ is given by: $w_{1-2}^{xx}=w_{3-4}^{yy}=1/3, w_{3-4}^{xx}-w_{1-2}^{yy}=-1/6, w_{5-8}^{xx}=w_{5-8}^{yy}=-1/24, w_{1-4}^{xy}=0, w_{5-6}^{xy}=1/4, w_{7-8}^{xy}=-1/4$. This choice can reduce the spurious velocities at the interfaces \cite{pooley_2008}.

\subsection{The Colour Gradient lattice Boltzmann Model}
An immiscible fluid model developed from Lattice Gas Cellular
Automata was introduced by Gunstensen and Rothman \cite{color1991}. The
particles in this model are coloured either red or blue and therefore it
is normally called the Colour Gradient Method. The surface tension is
introduced by adding a perturbation to the collision operator while
 adhering to the Navier-Stokes equations in homogeneous
regions. A recolouring step is used after the surface tension
perturbation calculation in order to achieve zero diffusivity of one
colour into the other.\\

We use $f_i^r$, $f_i^b$ and $f_i$ to denote the distribution functions
of a red fluid, a blue fluid and their combination respectively. A perturbation is introduced to generate the surface tension. The
surface tension can be expressed as a local anisotropy in the
pressure: the pressure measured normal to the surface is larger than
that tangential to the surface. The pressure in a the LBM simulation is proportional to the density, so that surface tension can be generated by preferentially placing  particles in
directions normal to the interface rather than tangential. Mass and
momentum are conserved. The color gradient G is defined as:

\begin{equation}
G_c(x,t)=\sum\limits_i e_i \sum\limits_j(f_j^r(x+e_i\Delta
t,t)-f_j^b(x+e_i\Delta t,t))
\end{equation}

The perturbation of the populations is given by:
\begin{equation}
f''_i(x,t) = f'_i(x,t) +A|G_c(x,t)|(\frac{(e_i\cdot G_c)^2}{G_c \cdot G_c}-\frac{1}{2})
\end{equation}

where the parameter $A$ determines the strength of the perturbation. On the other hand it also determines the strength of the surface tension. The surface tension $\sigma$ is given by $\sigma=\frac{192\rho A}{\omega}$, where $\omega$ is the inverse of the relaxation time \cite{rothman_1997}.\\


A redistribution of colour forces the particles to move towards the
regions occupied by particles with the same colour. This recolouring step
enables us to achieve separation of the fluids. It is carried out
by the following maximization problem \cite{color1991}:

\begin{equation}
W(f_i^{r''},f_i^{b''})=\max\limits_{f_i^{r''},f_i^{b''}}[\sum\limits_i
(f''_{r_i}-f''_{b_i})e_i]
\end{equation}

subject to the following constraints:

\begin{eqnarray}
&\rho''_r=\sum\limits_i f_i^{r''}=\rho_r&\\
&f_i^{r''}+f_i^{b''}=f_i&
\end{eqnarray}

A general colour gradient lattice Boltzmann model can be summarised
as:

\begin{itemize}
\item Single phase collision.
    %\begin{equation}
    %f_i=f_i^r+f_i^b
    %\end{equation}
\item Add a surface tension perturbation to $f'_i$ obtaining $f''_i$
\item Recolouring
\item Streaming
\end{itemize}


\section{Verifications for the multi-component lattice Boltzmann method}
\subsection{Phase separation simulation}
With the multi-phase lattice Boltzmann models, we can simulate phase separation and its dynamics. A $200 \times 200$ domain with an average density of $200 mu lu^{-1}$ ($mu$ is a lattice mass unit, and $lu$ is a lattice length unit) is computed with the Shan-Chen multi-phase lattice Boltzmann model (Figure \ref{PhaseSeparation}). The initial density is specified with random variations to prevent a metastable situation. The surface area is minimized as a consequence of the free energy minimization. This validation test is also simulated using the Free Energy LB model and the Color Gradient model, for both models, the mixture separated into two immiscible phases, which is similar to Figure (\ref{PhaseSeparation}).  %In the liquid-vapor system, condensation and evaporation phenomena can also be discovered during the simulation. This is due to the relatively high vapour density.
\begin{figure}[H]
\begin{center}
\includegraphics[width=5in]{phaseseperation.eps}
\end{center}
\caption{Liquid-vapor phase separation dynamics at time 0, 100, 200, 400, 2000, 6000 $ts$}\label{PhaseSeparation}
\end{figure}

\subsection{Simulation of wetting surface}
Boundaries with different wettabilities can be simulated with multi-component LBM. Droplets on a wetting surface with contact angles $20^\circ,90^\circ,150^\circ$ are simulated by the Free Energy model and are shown in Figure (\ref{wettingb}). The value of the contact angle can be set in advance using Equation (\ref{contact_angle_equation},\ref{wettingFE}). The Shan-Chen model and the Color Gradient model can also simulate binary fluid with wetting boundary and the results are in good agreement with analytical solutions.

 \begin{figure}[H]
\begin{center}
\includegraphics[width=3in]{wetting.eps}
\end{center}
\caption{Droplets simulated on wetting surface with contact angle $20^\circ$(left),  $90^\circ$(middle), $150^\circ$ (right). The liquid is represented in red while the vapour phase is shown in blue.}\label{wettingb}
\end{figure}

\subsection{Poiseuille flow simulation for an immiscible binary fluid system with viscosity contrast}
2D Simulations of two fluids with a viscosity ratio of 100 in a channel
have been carried out and compared with theoretical predictions.(Figure\ref{CGFE_v}) The densities of the two components are kept the same. The system is initialized with
substance 0 in the middle and substance 1 on both sides near the
walls. An initial density value of 1 is applied to substance 0 and 1, the initial viscosity of two components were set as 0.2 and 0.002 in order to achieve a viscosity ratio of 100. \\

\begin{figure}[H]
\begin{minipage}[t]{0.45\linewidth}
\centering
\includegraphics[width=\textwidth]{CG_poisseuille.eps}
%\caption{Simulated velocity of Color Gradient Model \label{CG_v}}
\end{minipage}
\hfill
\begin{minipage}[t]{0.45\linewidth}
\centering
\includegraphics[width=\textwidth]{FE_poisseuille.eps}
%\caption{Simulated velocity of Free Energy Model \label{FE_v}}
\end{minipage}
\caption{Simulated velocity of the Color Gradient Model (left) and the Free Energy Model (right) with viscosity ratio 100} \label{CGFE_v}
\end{figure}

The simulation results of both the Free Energy Model and the Colour Gradient Model 
give excellent agreement with the analytical solution. We should note that the interface thickness of the Free Energy Model is around $6$ lattice
units and the interface location is found at $x=38$ instead of the
initial position of $x=35$. An interface movement of 3 lattice units is discovered from the results. The interface thickness of the Color Gradient Model is around 4, which is smaller than that of the Free Energy Model. An interface shifting is observed in the colour gradient model of
around 1.5 lattice units which is significantly smaller than the
free energy model. The recolouring algorithm separates the two immiscible fluids very well
and gives nearly 0 diffusivity as expected. The Shan-Chen model failed to predict the velocity, a unfavoured diffusion near the interface of two components was observed. The big diffusion can be eliminated by using big inter-particle force, however, this will decrease the numerical stability significantly. In our numerical experiments, we always found diffusive interfaces even for very big inter-particle forces which caused the simulations very unstable.

   

\subsection{Capillary fingering simulation}
Capillary fingering is a kind of hydrodynamic instability which appears in various displacements during oil/gas production. One fluid is displaced by another fluid with lower viscosity, along a
channel with non-slip walls. A growing finger of the driving fluid will be produced if the capillary
number $Ca$ is big enough.


\begin{equation}
Ca=\frac{u_t \rho \nu_2}{\sigma}
\end{equation}

where $u_t$ is the velocity of the tip of the finger, $\nu_2$ is the
viscosity of the driving fluid and $\sigma$ is the surface tension.


Fingering in a channel has been investigated using two different lattice Boltzmann binary models: the Free Energy Model and the Color Gradient Model. Figure (\ref{fingering2}) shows the evolution of fingers simulated by the Free Energy Model. From top to bottom, we show the finger evolution for surface tensions of 0.0678, 0.0389, 0.01985, 0.00992 in lattice units respectively. No finger will be produced if the surface tension is high. When the surface tension decreases, fingers are observed. Halpern and Gaver \cite{halpern_1994} studied the fingering phenomenon in a channel.  The width of the fingers produced by different capillary numbers $Ca$ are measured. An empirical formula was proposed to describe this relationship between the width of capillary fingering and capillary number. Their results are plotted in Figure (\ref{fingering1}) with a solid line. A relative width value for the width of the channel was used. Our results from the Free Energy Model and the Color Gradient Model are shown in the same figure using triangle points and star points respectively. We only made one data point for high capillary number because the simulation became unstable as the capillary number increase due to the increasing velocity. On the other hand, the data of high capillary number is not of interest for reservoirs engineering, because the capillary number of flow in the reservoir is normally very small.\\

\begin{figure}[H]
\begin{center}
\includegraphics[width=2.5in]{Capillary_Fingering1.eps}
\end{center}
\caption{Fingering evolution with different surface tension at a time interval of 1000 time steps. From top to bottom, we show the finger evolution for surface tensions of 0.0678, 0.0389, 0.01985, 0.00992 in lattice units respectively.}\label{fingering2}
\end{figure}



\begin{figure}[H]
\begin{center}
\includegraphics[width=4in]{fingering.eps}
\end{center}
\caption{Fingering as a function of Capillary number. The results from the Free Energy and Color Gradient Model are represented with triangle and star points, respectively}\label{fingering1}
\end{figure}

Good agreement is achieved although some small discrepancies are found. These discrepancies might come from the boundary conditions. These numerical examples show that both the Free Energy Model and the Color Gradient Model are capable of simulating the capillary fingering phenomenon. It is worth mentioning that the maximum viscosity ratio for dynamic interfaces simulations is around 20. A higher viscosity ratio may lead to numerical instability and sometimes crashes the code.% Knowing this value can help people to identify if the Free Energy Model or the Color Gradient Model is suitable for their specific case.

\section{Parallel LBM implementation}
The LBM simulator is able to simulate various complex flows with extremely complicated boundaries. However, this is very time-consuming even with modern computers. 
%Recent development of computer-aided imaging enable us to obtain high resolution 3D models of rocks. The porosity of these voxel models are normally less than $40\%$, 


In a single phase flow simulation, most of the computation is local. In the collision step, the distribution of fluid particles will be changed and no distribution of neighbouring points is needed. The situation is the same for the calculation of macroscopic quantities. Communication between points is only needed in the streaming step at which the updated distribution functions will stream to the adjacent mesh points. The direction of streaming is determined by the velocity of the particles on the mesh points. Parallel computing has been available for decades and could be a solution to improve the efficiency of the LBM code.\\

Several studies have been carried out to investigate the efficiency of the parallel lattice Boltzmann method. Skordos (1995) compared two parallel CFD methods, finite differential method and lattice Boltzmann method. An equal partition strategy was used in this study. He investigated the relationship between problem sizes, size of each partition, number of partitions and the performance of the code. A $50\%$ increase of efficiency is obtained if the number of processors is more than 15. Martys et al. (1999) implemented a parallel lattice Boltzmann code to simulate multiphase flows. A speedup slightly lower than linear was achieved. \\

In this work, a Message Passing Interface (MPI) library was used to carry out the communications between partitions. The computational domain is divided equally into several partitions in the x direction. This is shown in Figure (\ref{partition1}). Particles from different partitions are illustrated with different colours. The particles in the dotted line area will stream to the adjacent nodes and therefore need to communicate with particles that are stored on another processor. An MPI subroutine has been developed to exchange distribution function values for particles in this area.     

\begin{figure}[H]
\begin{center}
\includegraphics[width=3in]{partition.eps}
\end{center}
\caption{Schematic diagram of the partition geometry}\label{partition1}
\end{figure}



\section{Summary}
The work that has been accomplished can be summarised as:

\begin{itemize}
\item Implementation of a 2D and 3D parallel lattice Boltzmann method code for single phase flow simulation. The code is capable of:
	\begin{itemize}
	\item Simulating single phase fluid flow with Reynolds numbers of up to 1000
	\item Simulating single phase fluid flow within a complex geometry without viscous boundary effects by using the MRT scheme
	\item Computing absolute permeability of the porous medium
	\end{itemize} 
\item Implementation of a 2D and 3D lattice Boltzmann method code for an immiscible binary fluid system or a single-component multi-phase system. The code is capable of :
	\begin{itemize}
	\item Simulating multi-phase flow with a density ratio of up to 20 in porous media using the Shan-Chen pseudo potential model
	\item Simulating an immiscible binary fluid system with a viscosity ratio up to 20 in porous media using the Free-Energy model and the Colour Gradient model 
	\item Simulating wetting effects of the surface
	\item Simulating interfacial dynamics of an immiscible binary fluid system 
	\end{itemize}
\item Verification work of the code package
\item Absolute permeability calculation for pipe flow and flow in fibrous porous media
\end{itemize}

The original work I have done is:
\begin{itemize}
\item Comparison of three multi-component lattice Boltzmann models with numerical experiments; the strengths and limitations of these models were discussed
%\item Absolute permeability calculation for pipe flow and flow in fibrous porous media using the  lattice Boltzmann method 
\end{itemize}

Based on our preliminary results, we believe that the lattice Boltzmann method is a promising tool for the estimation of transport properties of fluids in porous media due to the following reasons:
\begin{itemize}
\item Chapman-Enskog expansion analysis showed that the lattice Boltzmann method can recover the Navier-Stokes equation at macroscopic scale\cite{ekexpansion}.
\item Young-Laplace law behavior is reproduced correctly.
\item Numerical experiments have shown that wetting phenomena can be simulated properly.
\item Poiseuille flow simulation showed that the lattice Boltzmann method can estimate the velocity very accurately even if only one fluid node is imposed between solid surfaces.
\item Poiseuille flow simulation for a binary fluid system with viscosity contrast showed that the kinetic behaviour can be simulated accurately with a tiny interface shifting. 
\item Simulation of flow in fibrous porous media showed that the LB method can estimate the absolute permeability accurately within a complex porous geometry.
%\item Absolute permeability simulation and the experiments showed that the lattice Boltzmann method can predict the absolute permeability of the rocks with reasonable accuracy if a fine enough mesh is imposed in the simulation.  
\end{itemize}



\chapter{Plan for future work}
\section{General plan for future work}
Initially, we will continue to  improve the efficiency of our code. The reliability of our parallel single-phase code has been verified by our numerical experiments. From a preliminary study of the absolute permeability of two sandstones (Clashach and Doddington), we learned that a very fine mesh is required to achieve reasonable accuracy. A standard storage strategy is used in the code package at the moment. The data on all the nodes are stored and the computing operation is carried out on all these nodes. Many rocks of interest have a porosity ranging from $8\%$ to $40\%$. Only these $8\%$ to $40\%$ fluid nodes take part in the computing operations in the lattice Boltzmann method. The standard storage strategy wastes a large amount of memory and computing time for solid nodes which are useless in the simulation. We will introduce a sparse storage strategy based on Mattila et al (2007)\cite{mattila_2007} into our parallel lattice Boltzmann method code package. In this sparse storage strategy, only the fluid nodes will be stored which can save a significant amount of memory and computing time. The adjacent node information including coordinates of adjacent nodes and physical address of adjacent nodes will be compressed and stored properly for the communication operations in the streaming step. According to literature data, we expect to achieve a reduction of $30\%$ to $76\%$ for memory usage along with an improvement of efficiency by a factor of 3-4. \\

Several numerical tests including Poiseuille flow simulation, the simulation of flow past a cylinder, the calculation of absolute permeability in fibrous porous media and calculation of flow in real rock geometries will be carried out to verify the reliability of this sparse storage strategy. The performance will then be compared with that of the lattice Boltzmann method using the standard storage strategy. When the validity of this sparse storage strategy has been confirmed, it will also be applied to our multi-component LBM code. \\

We will parallelize our 3D multi-component LB code to improve the efficiency. A partition strategy similar to the one used in our single-phase lattice Boltzmann code will be imposed.  A parallel implementation along with sparse storage strategy enables us to simulate geometries with large size. We expect to be able to simulate a porous medium with a grid size of 1500x1500x1500. Some spurious currents were found in our multi-component lattice Boltzmann simulations. These spurious currents in the Free Energy Model can be due to an inappropriate choice of distribution functions. Pooley et al (2008)\cite{pooley_2008} proposed a modified distribution function for the Free Energy Model so that the spurious velocities may be decreased by an order of magnitude compared to previous models. We are going to implement this new scheme in our Free Energy lattice Boltzmann code for multi-component flow simulations to improve the accuracy of the code. For the Shan-Chen model, these spurious velocities can be due to the insufficient isotropy of the discrete gradient operator\cite{shan2006}. We are going to implement finite difference gradient operators with a higher order of isotropy in the multi-phase and multi-component Shan-Chen LBM simulator to reduce the spurious currents.\\

Secondly, we are going to verify our single-phase and multi-phase multi-component LB code with experimental data. Several plastic models with capillary channels and junctions with different geometries will be made to study the capillary trapping of $CO_2$ and the effect of surface wettability on the trapping process, both experimentally (by E.Chapman) and numerically (this PhD project). Two example geometries have been given in Figure(\ref{Junction}) in Chapter(\ref{chap1}). Experiments of fluid flow in these plastic models will be carried out in a related separate PhD project (E. Chapman). The corresponding numerical simulation will be carried out with our code. The simulation results will be directly compared with the experimental results. \\

A micromodel geometry based on a 2D thin section of a Sucrosic Dolomite was designed to study the transport phenomena, capillary trapping of supercritical $CO_2$ and the effect of surface wettability  (with E. Chapman). A digital model was established first to tune the porosity and permeability of the model to match the experimental values \cite{knackstedt_2007} using our single-phase lattice Boltzmann code. Following the computational design method, etched glass micromodels are now being produced. We will use the experimental models to measure the  absolute permeability and other transport properties. To our knowledge, this is the first time that a micro-fluidic micromodel has been designed by computation, prior to the real geometry being produced. We will verify our design and numerical results with the experiments.\\


Images are to be taken from several sandstone and carbonate rock samples using multi-scale imaging techniques including confocal microscopy and x-ray micro-tomography (a separate PhD project conducted by S.Shah). The 3D geometry of the rocks will be established with these images. These 3D geometries will be used as input files for our lattice Boltzmann simulator. Single-phase flow will be simulated to compute the absolute permeability. The Multi-Component lattice Boltzmann code is to be used to mimic a primary drainage displacement and an imbibition process by injecting non-wetting fluid into the sample which is saturated with wetting fluid.  Periodic boundary conditions will be applied on the inlet and outlet. The wetting fluid exiting the geometry will re-enter at the inlet. The relative permeability will be calculated and compared with experimental data. We will simulate the primary drainage and imbibition of oil and brine or oil and supercritical $CO_2$ in the rocks. In addition, we will also study the steady state flow of binary fluid systems in rocks. To simulate steady-state flow, wetting fluid and non-wetting fluid are randomly distributed in the geometry to mimic a target saturation value. When the flow reaches steady state (the change rate of the velocity field is less than $10^{-6}$), the relative permeability will be computed. Experimental data will be used again to evaluate the performance of our code for estimating transport properties within reservoir rocks.\\

Finally, to simulate the flow and transport phenomena on a core scale, we will implement a lattice Boltzmann method coupled with Darcy's law \cite{dardis_1998}. In this model, the porous medium will be parameterized in terms of  spatially varying porosity and permeability. The collision algorithm will be modified so that the flow is governed by Navier-Stokes equations and Darcy'Law simultaneously. This core scale simulation enables us to estimate the distribution of water, oil and supercritical $CO_2$ at the reservoir scale. Such data will help the operational design of reinjecting produced CO2.\\

\section{Summary of future work and corresponding actions}
\begin{itemize}
\item Implementation of sparse storage data structure.
	\begin{itemize}
	\item A sparse storage data structure is to be developed based on the algorithm proposed by Mattila et al. \cite{mattila_2007}
	\item The code will be tested with several numerical experiments. The performance will be compared with the LBM code with standard storage strategy\cite{shan2006}.
	 
	\end{itemize}
\item Parallel implementation of multi-phase multi-component LBM code
\item Reduction of spurious currents in multi-component LBM simulation
	\begin{itemize}
	\item The spurious currents in the Free Energy model are expected to be reduced by using a novel distribution function proposed by Pooley\cite{pooley_2008}
	\item The spurious currents in the Shan-Chen model are to be reduced by implementing finite difference gradient operators with a higher order of isotropy
	\end{itemize}


\item Verification of the code. Several verification tests will be carried out to verify the reliability of the code. Theoretical predictions are to be used to evaluate the performance of the code. The tests will include:
	\begin{itemize}
	\item Poiseuille flow in a channel
	\item Karman street vortex shedding simulation
	\item Flow past a cylinder
	\item Flow past arrays of spheres
	\item Flow in fibrous porous media
	\item Young-Laplace law verification (By simulating a series of bubbles with different radius, the pressure difference will be measured to verify the Young-Laplace Law)
	\item Droplets on surfaces with different wetability
	\item Capillary rise of an immiscible binary fluid system in a capillary tube
	\item Capillary fingering  of immiscible binary fluids with viscosity contrast
	\item Poiseuille flow simulation of an immiscible binary fluid system with viscosity contrast in a channel
	\item Snap-off phenomena simulation
	\end{itemize}

\item The lattice Boltzmann method code is to be used to simulate micro-fluidic micro-models with capillary channels and junctions. Flows and effect of wettability of surface will be studied in order to study the capillary trapping of $CO_2$ and the effect of surface wettability on the trapping process.

\item The flow in a porous medium geometry based on Sucrosic Dolemite is to be simulated. Experimental data will be used to evaluate the performance of the code.

\item Flow in Sandstone and carbonate rocks will be calculated. The geometry is to be taken using multi-scale imaging techniques. The dynamic displacement process including a primary drainage displacement and an imbibition process will be studied using our lattice Boltzmann code. Absolute and relative permeability is to be computed and compared with experimental data.

\item A lattice Boltzmann model coupled with Darcy's law will be developed to simulate the flow in the reservoir at core scale. The algorithm for this implementation will be based on the model proposed by Dardis and McCloskey \cite{dardis_1998}

\section{Time Plan}
\begin{figure}[H]
\begin{center}
\includegraphics[width=5in]{timetable.eps}
\end{center}
\caption{The time plan}
\end{figure}


\end{itemize}




\bibliographystyle{ieeetr}
\bibliography{ref}



\appendix
\input{appendix1}
\input{appendix2}

\end{document}
